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ABSTRACT 

We present an improved formula for the anomalous radial velocity of the star dur- 
ing planetary transits due to the Rossiter-McLaughlin (RM) effect. The improvement 
comes from a more realistic description of the stellar absorption line profiles, taking 
into account stellar rotation, macroturbulence, thermal broadening, pressure broaden- 
ing, and instrumental broadening. Although the formula is derived for the case in which 
radial velocities are measured by cross-correlation, we show through numerical simu- 
lations that the formula accurately describes the cases where the radial velocities are 
measured with the iodine absorption-cell technique. The formula relies on prior knowl- 
edge of the parameters describing macroturbulence, instrumental broadening and other 
broadening mechanisms, but even 30% errors in those parameters do not significantly 
change the results in typical circumstances. We show that the new analytic formula 
agrees with previous ones that had been computed on a case-by-case basis via numerical 
simulations. Finally, as one application of the new formula, we reassess the impact of 
the differential rotation on the RM velocity anomaly. We show that differential rotation 
of a rapidly rotating star may have a significant impact on future RM observations. 
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Introduction 



Transiting exoplanetary systems provide valuable opportunities to learn about the nature 
of the exoplanets, their orbits, and their host stars. In particular, when we measure the radial 
velocity (RV) of a star during a planetary transit, we see an anomalous Doppler sh ift (beyond the 
usual orbital RV) which is called the Rossiter-McLaug hlin (hereafter, RM) effect (lRossiterlll924l ; 
McLaughlinlll924l : lHosokawalll953l : lAlbrecht et al.ll2007l ). It arises because a portion of the rotating 
stellar disk is blocked by the planet. The partial occultation brings about a distortion in the 
spectral lines, which is manifested as an anomalous RV depending on the position of the planet on 



the stellar disk (see, e.gjQueloz et alJl20(Kl : lOhta et al.ll2005l : IWinn et al.ll2005l : iNarita et alJl2007l : 



Collier Cameron et al.ll2010l ). The time variation of the RV anomaly reveals the (sky-projected) 



angle A between the planetary orbital axis and the stellar spin axis. Measurements of this angle 
have proved to be an important observational clue to the origin of close-in giant exoplanets. 

It is widely assumed that close-in gas giants, of which more than 1 00 are known, formed at a few 
AU a way from their host stars and subsequently "migrated" inward (ILin et al.lll996l : iLubow Ida 
2010) . Many planetary migration scenarios have been investigated, and som e of them predict small 
values o f A ~ 0° while others allow larger spin-orbit misalignment an gles (e.g.. lFabrvckv fc Tremaine 



20071 : IWu et al J 120071 : iNagasawa et alJ 12008 : IChatteriee et alJ 12008 ). The observed distribution of 
A and its dependence on the host star properties (such as masses and age s ) can be important clues 
to understand the origin of close-in giant exoplanets (jWinn et al.l l2010al : iFabrvckv fc Winnl 120091 : 
Morton Johnsonlboill ). 



Observations of the RM effect have now become almost routine (IWinn 1 120101 : iMoutou et al 



20 111 ). However, it is important to remember that the relationship between the observed RV 
anomaly, and the position of the planet on the stellar disk, is not completely straightforward. This 
is because the RM effect is actually a spectral distortion, even though it is frequently studied 
as though it were a pure Doppler shift|j Many alternatives have been pursued to calibrate the 
relationship betw een t he observed sig nal and the underlying parameters of the planet and star. 
Ohta et al.l ( 20051 ) and iGimenea ( 200a ) derived analytic formulas for the RV anomaly based on the 
computation of the first moment (the intensity-weighted mean wavelength) of distorted spectral 
lines. This approach is simple and convenient because the compu ted velocity anomaly does not 
depend on the intrinsic shape of spectral lines (jHirano et al.ll2010l ). and has been useful for quick 
computations where high accuracy is not essential, and for gaining insight into the parameter 



dependence of the RM veloci ty anomaly. I Winn et al.1 (120051 ). however, noted that the analytic 
formula by lOhta et al.l (|2005l ) (hereafter, the OTS formula) deviates from the results based on 
a more realistic numerical calibration; they simulated spectra exhibiting the RM effect for many 
different positions and sizes of the planet, and then analyzed the mock spectra with the same 
data analysis codes that are used routinely to derive precise RVs with the High Resolution Echelle 



An alternative is to m odel the line profiles directly, as has been done by lAlbrecht et al.l ( 20071 ) and 
Collier Cameron et al. (2010). which can be advantageous in some circumstances. 
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Spectrometer (HIRES) installed on the Keck I telescope. As a result, they showed that the OTS 
formula disagrees with the numerical calibration by about 10% in terms of the RM amplitude for 



the ca se of HD 209458. For this reason, subsequent studies (e.g., IWinn et alj|2005l : iNarita et al 



2009al ) have relied upon numerical calibration of the relation between the anomalous RV and the 



position and size of the planet, which is done on a case- by-case basis depending on the stellar 
parameters. 

It would be more convenient to rely on a single analytic formula than to perform these laborious 
numerical simulations for each syste m. An analy t ic tre atme nt also proyides in sight into the reason 
for the limitation of the formulas of lOhta et al.l (|2005l ) and iGimenea (|2006l ) . iHirano et al .1 (|2010l ) 
took a step in this direction, pointing out that the discrepancy between the OTS formula and the 
simulated results was a consequence of the algorithm used to estimate the RV anomaly. While 
the OTS formula was derived by computing the first moment of the distorted line, in practice 
the RVs are computed by cross-correlating an observed spectrum with a template sp ectrum of the 



same star [e.g., for the High Accuracy Radial velocity Planet Searcher (HARPS); iTriaud et al 



(|2009l )]. or by forward-modeling to fit an observed spectrum with superimposed iodine absorption 
lines [e.g., f or the High Dispersi o n Spe ctrograph (HDS) on the Subaru telescope and Keck/HIRES; 



Sato et al.l (|2002l ); iButler et all (1996T)] . Using a simplified description of a single spectral line (a 



Gaussian function) , IHirano et al .1 (|2010l ) compared the RV anomalies derived by computing the first 
moment and by cross-correlation. They showed that the two methodologies yield different velocity 
anomalies, in a manner that qualitatively explains the previous numerical findings. In particular 
they showed that the deviations between the OTS formula and the results of cross-correlation are 
larger for more rapidly rotating stars. 

In this work, we take the next step by developing a more realistic description of stellar line pro- 
files, in order to derive a more accurate analytic formula. Instead of using a simple Gaussian model 
for a spectral line profile, we include realistic kernels for rotational broadening, macroturbulence 
and other effects such as instrumental broadening due to the finite resolution of a spectrograph. 
We test and validate the new analytic formula through various numerical simulations, and show 
that it is accurate enough for the real data analysis. 

This paper is organized as follows. In Section [2j we derive the new analytic formula for the 
RM effect assuming an analytic function for the stellar line profile. We present the definitions 
and the result there, while the detailed derivation of the main finding is described in Appendices 
lAl and [Bj In order to make sure that the new analytic formula is a good approximation for the 
observed velocity anomaly due to the RM effect, we compare it with numerical simulations using 
mock transit spectra in Section [3l Also, we check on the magnitude of systematic errors due to 
imperfect knowledge of the parameters describing the absorption line profiles for a given star. As 
an application of the new analytic formula, we try to reassess the impact of stellar differential 
rotations on the RM velocity anomaly taking the XO-3 system as a test case in Section [H The 
final section ($5]) is devoted to discussion and summary. 
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2. Derivation of the New Analytic Formula for the RM Effect 



In this section, we derive the new analytic formula that describes the velocity anomaly during 
a t ransit. We begin wi th our description of the stellar absorption lines. We follow the formulation 
bv lHirano et al.l (|2010l ) but slightly change the basic equations in order to describe the stellar line 
profiles more realistically. Since the velocity field on the stellar surface is of primary importance, it 
is more convenient to express all the functions in terms of velocity rather than wavelength. In what 
follows, the velocity component v indicates the velocity shift relative to the center of an absorption 
line. This is related to the wavelength shift A A by the usual formula AA/Ao = v/c, where Ao is 
the central wave length of th e absorption line and c is the speed of light. Following the model of 
spectral lines by [Gray (|2005l ). we write a stellar line shape J r s tar(v) as 



T sUr {v) = -S(v)*M(v), 



(1) 



where S(v ) is the intrinsic stellar line shape in the absence of stellar rotation and macroturbulence 
(for which we will give an explicit expression later), and M(v) is the broadening kernel due to 
stellar rotation and macroturbulencdE The symbol * indicates a convolution between two functions. 
Since the continuum level and the normalization factor in the spectrum do not affect the result in 
estimating the velocity anomaly during a transit, for convenience we subtract the continuum level 
so that Fst&riv) becomes zero in the limit of v — > ±oo. Furthermore we normalize the spectrum so 
that 



Fstnr^dv 



-1. 



(2) 



The minus sign indicates that J- S ta.r(v) describes an absorption line. The rotational-macroturbulence 
broadening kernel M{v) is calculated by disk-integrating the Doppler-shift component of the stellar 
surface due to both stellar rotation and macroturbulence. We adopt "the radial-tangential model" 
for macroturbulence, for which the kernel in the absence of rotation is 



Q(v) 



1 



2^ 



1 



£ cos t 



+ 



1 



, C sin 6 t 



(3) 



C cos 6 C sin ( 

where £ is the macroturbulent velocity param eter and 6 is the angle between our line-of-sight and 
the normal vector to the local stellar surface (|Gravll2005l . page 433). The angle 6 is related to the 
coordinate (x, y) on the stellar disk by 



cos I 



x 2 + y 2 



R 2 



sin I 



\Jx 2 + y 2 



(4) 



where the y-axis is taken to be along the sky projection of the stellar spin axis, and R s is the radius 
of the star. Assuming a quadratic limb-darkening law, the disk-integrated line broadening function 



2 We here assume a symmetric line profile and ignore the convective blueshift (CB) effect, discussed by 
Shporer fc Brownl (|201lh . 
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due to stellar rotation and macroturbulence is expressed as 

1 - u\(l - cos#) — 1*2(1 — cos#) 2 



M(v) 



entire disk 



7r(l - ui/3 - 1*2/6) 



Q(v — xQsinig) 



dx dy 



(5) 



where u\ and 112 are the limb-darkening coefficients, $7 is the angular spin velocity of the star, and i s 
is the inclination angle of the stellar spin axis relative to the line of sight (lGrayll2005l ) . The Doppler 
shift — x^sin i.s in t h e fun ction Q(v) is the consequence of stellar rotation, neglecting differential 
rotation. As iGrayl (120051 ) pointed out, the broadening kernel M(v) cannot be expressed as a 
convolution of the two different broadening kernels of the stellar rotation and the macroturbulence. 
As we will show, the coupling between rotational broadening and macroturbulent broadening plays 
an important role in estimating the velocity anomaly due to the RM effect, especially when the 
macroturbulent velocity is appreciable when compared to the rotational velocity of the star (see the 
difference between line profiles with and without macroturbulence shown in Figure [1]). Indeed, this 
co upling betwe e n rota tion and macroturbulence was neglected in the previous numerical calibrations 
by I Winn et all (|2005l ) and others. 
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Fig. 1. — A schematic plot of the line profile during a planetary transit. Each line has a different 
broadening kernel. For visual clarity, the line profiles are vertically shifted by 0.2 along the intensity 
axis. Black line: an intrinsic line profile without stellar rotation and macroturbulence, described 
as a single Gaussian function with standard deviation f3 = 1 km s" 1 . Red line: after convolving 
with a pure-rotational broadening kernel (no macroturbulence), with vs'mi s = 6 km s _1 . Blue 
line: after convolving with a rotational-macroturbulent broadening kernel with t;sini s = 6 km s _1 
and £ = 4 km s _1 . For the latter two cases (red and blue lines) the spectral contribution from 
the portion occulted by the planet has been subtracted from the profiles, assuming a planet with 
(R p /R s ) 2 = 0.01. The line profile with macroturbulence (blue) has elongated wings and the transit 
signal is nearly invisible. 



Next, we compute the line shape during a planetary transit. During a transit, the spectral 
contribution of the portion blocked by the planet is written as 

F p i anct (v) = -S(v)*M'(v), (6) 

where M'(v) indicates a kernel similar to that given in Equation ([5]) but for which the disk- 
integration should only be performed over the blocked part of the stellar surface, rather than the 
entire stellar disk. As long as the planet is sufficiently small relative to the star, the Doppler shift 
— xf2sini s in Equation ([5]) is nearly constant over the integration region. Thus, if we define X as 
the x-coordinate of the intensity-weighted center of the eclipsed portion of the star, we can remove 



-7- 



the macroturbulence kernel 0(t>) from the integral and define the following two useful quantities: 
j- ff 1 — ni(l — cos 6*) — U2(l — cos 6*) 2 dx dy 

J J occulted portion 7t(1 - Ul/3 - 7X 2 /6) R 2 S 

v p = XQs'mig, (8) 
so that J^ianet (v) becomes 

Fpl a nct(v) = -f S(v)*e(v-V p ). (9) 

The first quantity, /, is the instantaneous fractional decrease in flux due to the transit. The second 
quantity, v p , is the rotational radial velocity of the occulted portion of the stellar disk, which is 
occasionally referred to as the "subplanet velocity." With these definitions the stellar line profile 
during a transit is expressed as 

^transit (") = -^star(«) - -^planctO) = S(v) * M(v) + fS(v) * Q(v - V p ). (10) 

It should be noted that the macroturbulent kernel @(v) remains in the modeled transit line profile. 
This treatment is necessary since the two effects of rotational broadening and macroturbulence are 
coupled with each other. In short, the lin e profile during a tr ansit expressed by Equation (llOj) is 



different from the line profile modeled by iHirano et al.l (|2010l ) (Eq.[ll]) in two senses: Equation 



(HOP explicitly involves the effect of macroturbulence, and it is expressed in terms of velocity. 
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Table 1: Summary of symbols used in this paper. 



Symbol 


Meaning 


Typical Range 


f 
J 


tVip itmf flTitflnpniK frfipi"innfil HppTpflQp in flnv Hup 1~n 1~Vip frfnr^il - iT^n 171 1 

LillC lllol Ctll LcllH_.WLlo 1 1 CI Lill W 1 1 CI 1 L- 1 \_.(XO\_^ 111 11 L1.A. LI 11' L w lj 11*1. Ijldllollj \ LlivJ ■ III / 


00 - 02 


v~. 

Up 


the surinlanpf, vploritv (TOn \H\ ) 

uii*j o ui is yj lcaii.±.kj \j v V;iw v^i u y \ i lj i / 


-hi j sin ?,„ 

— i — yj oiii V § 


Av 


the velocity anomaly due to the RJVI effect 






tli6 limb-dcLrkGning pcLrcimGtGrs for tlic cjUcLclrcitic limb-dcLrkGning l&w 


0.3 - 0.5 


v sin i s 


the stellar spin velocity 




- 1 ert 


the stellar effective tenrnerature 




7 « 


th6 inclination of tli6 stellar spin axis measured from our line-of-sight 


0° - 90° 


/ 


the latitude on the stellar surface 


±90° 


L S 


the stellar radius 






the position of the transiting planet on the stellar disk 


±R S 


a 


the coefficient of differendtial rotation 


±0.02 


p 


the Gaussian dispersion of spectral lines (Eq. j20j) 


2.5 - 4.5 km s" 1 


7 


the Lorentzian dispersion of spectral lines 


0.5 - 1.5 km s" 1 


c 


the macroturbulence dispersion 


2.0 - 6.5 km s" 1 


e 


the angle between the line-of-sight and the normal vector to the stellar surface 


0° - 90° 


A 


the spin-orbit misalignment angle 


±180° 




the microturbulence dispersion 


0.0 - 2.0 km s _1 




the frequency in the Fourier domain 




ft 


the angular velocity of the stellar spin 





Armed with the preceding results, we now express the velocity anomaly Av during a transit 
in terms of t he fra ctional flux decrease / and the subplanet velocity v p . Basically, we follow 



Hirano et al.l (|2010l ) in order to compute the best-fit value for the anomalous RVs; they cross- 



correlated the spectrum during a transit with a stellar template spectrum, and then calculated the 
best-fit value for the velocity anomaly Av by maximizing the cross-correlation function C{x): 

dC{x) 



dx 
C(x) 



x=Av 
oo 



= 0, 

Fstariv 



transit 



[v)dv. 



(11) 

(12) 



To proceed further, we need a specific model for the intrinsic line shape S{v). We here adopt the 
Voigt function for S(v): 



S(v) 



V(v;P, 7 ) 
G(v;f3) 

L(v;j) 



G(v;(3)*L(v; 7 ), 

1 _„,2/fl2 



-V 2 //3 2 



1 7 

7T V 2 + 7 2 



(13) 
(14) 

(15) 
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where (3 is the thermal velocity parameter and 7 is the Lorentzian velocity parameter (due to 
pressure broadening or natural broadening). These parameters are related to individual stellar 
properties such as the effective temperature, surface gravity, and the nature of each absorption 
line. Some line profiles of especially strong absorption lines (such as the Na D lines) are saturated 
and intrinsically different from the Voigt function in shape. However, most of the lines in the 
wavelength region used in RV analyses are relatively weak, by design, and are well approximated 
by the Voigt function in the absence of the stellar rotation and macroturbulence. 

Substituting Equations ([I]) and (fTU|) into Equations (fTTj) and (fT2j) . we compute the velocity 
anomaly Av due to the RM effect. Since further calculations are mathematically complicated, we 
describe the detail of the derivation in Appendix lAl and write down the result alone: 

POD 

/ exp(-2vr 2 /3 2 o- 2 - 4n7a)M(a)@(a) sm(2irav p )ada 

Av ^-4= roo J ° : : > (16) 

J exp(-2 7 r 2 ( S 2 a 2 - A^a)M(a) |m(<t) - f@(a) cos(2irav p )j a 2 da 

where 

/•oo 

M{v)e~ 27Tiav - 

00 



M{a) = I M{v)e- Zmav dv 

- 1 1 - u x (i - VT^fi) - u 2 (i - VT^fif 



1 — u\/3 — U2/6 

x { e — 'C 2 -^ 1 "* 2 ) + e -" 2c2CT2 * 2 }Jo(2vr^sin^t)t ( it, (17) 
0(<7) = - [exp{-«cos^) 2 a 2 } +exp{-(vrCsin0) 2 ( T 2 }] , (18) 



1 



where J n ( x ) is the Bessel function of the first kind (see also Appendix [B]) . Equation (|16p is the 
main finding in the present paper and will be used in the comparison with simulated results. 



3. Validity of the Analytic Formula 

3.1. Comparison with Numerical Simulations for Subaru/HDS 

One of the major differences between the derivation of the analytic formula (Eq. [16]) and the 
manner in which RV data are actually analyzed is that the analytic formula is based on the cross- 
correlation method while the data analysis (for Subaru/HDS and Keck/HIRES at least) adopts 
the forward modeling method using the Iodine cell for a precise wavelength calibration. Another 
difference between them is that the analytic formula assumes only one absorption line in its deriva- 
tion while the actual stellar spectra have many lines differing in depth, shape, and so on. Thus, in 
order to test the validity of Equation (fTBT) . we perform the mock data simulation described below 
and compare the results with the analytic formula. 
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3.1.1. Mock Data Simulation 



The mock data simulation with the Iodine RV calibration is described in detail by I Winn et al 



(|2005l ) and lNarita et al .1 (|2009al ). In order to generate mo ck spectra during a transit, we begin with 
theoretically synthesized spectra by ICoelho et al.l (|2005l ). whose spectral lines are modeled by in- 
corporating thermal broadening (including microturbulence) , and Lorentzian (natural or pressure) 
broadening. Since intrinsic line profiles depend on the stellar type, we first test the comparison for 
a GO star, whose effective temperature is 6000 K. When we obtain the synthetic spectrum from the 
synthetic spectrum library, we assume a representative value of the surface gravity (log g = 4.0) for 
a GO star with a planet and the solar abundance for metallicity ([Fe/H]=0.0). 

Using the synthetic spectrum, we generate the mock spectra during a transit by the following 
procedure. 



1. We broaden the synthetic spectrum by convolving the rotational-macroturbulent kernel M(v) 
(Eq. [5]), so that it represents the template spectr um of the a ctual star. We adopt the 
limb-darkening parameter s u-\ = 0.43 and ui = 0.31 (jClaretl I2004J ) . and the macroturbulence 
parameter £ = 4.3 km s _1 ( Valenti et al.ll2005l ). We try three different values for the rotational 
velocity of the star: usini s = 2.5 km s _1 , 5.0 km s , and 7.5 km s . 

2. We Doppler-shift the original unbroadened spectrum by v p , multiply by /, and then convolve 
the spectrum with the macroturbulence kernel Q(v) (Eq. [3]), so that the resultant spectrum 
represents the spectral contribution from the portion occulted by the transiting planet. For 
each of the three values of i>sini s , we consider / = 0.004, 0.008, 0.012, 0.016, and 0.020. 
We then assign 21 different values to v p evenly spaced from — vsini s to +t>sini s , yielding 
in total 105 different points in the (f,v p ) grid for each value of vs'mi s . The macroturbulent 
broadening depends not only on the x-coordinate but also on the y-coordinate on the stellar 
disk (see Eqs. [3] and [1]). For simplicity, however, we assume y = when we generate the 
mock transit spectra. 

3. We generate the mock transit spectra by subtracting the Doppler-shifted and scaled spectra 
in the second step from the broadened spectra created in the first step. 

4. We multiply the mock transit spectra by the iodine transmission spectrum used for cali- 
bration, and convolve the Star+12 spectrum with the representative instrumental profile of 
Subaru/HDS for the case of the slit width being 0.8", corresponding to the spectral resolution 
of R ~ 45000 § 



We then take these 105 mock spectra for each of three values of v sin i s , and use them as inputs 
to the RV analysis routine for Subaru/HDS. The RV analysis using the iodine cell is described in 



3 The choices for the instrumental profiles used in the mock simulation are explained later. 
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detail bv lSato et al.l (12002 ) (for Subaru/HDS) and Butler et all (|1996h (for Keck/HIRES). For each 
point of the (/, v p ) grid, the RV analysis routine outputs a velocity anomaly Av due to the RM 
effect. 



3.1.2. Results for a GO Star 

We compare the results (f,v p ,Av) based on the mock data simulation with the analytic ex- 
pression (fTUj) , Figure [2] shows the comparison between the simulated data points (f,v p ,Av) and 
the analytic formula (Eqn. [IB]). The three different panels are for the three different values of 
vsini s : (top) 2.5 km s _1 , (middle) 5.0 km s _1 , and (bottom) 7.5 km s _1 . For each of the five 
values of /, the data points indicated by color symbols show the simulated results; / = 0.004 in 
black, / = 0.008 in red, / = 0.012 in blue, / = 0.016 in purple, and / = 0.020 in green. 

In order to draw the analytic curves, we need to choose values of /3 and 7. Although they are 
related to intrinsic stellar line profiles, each spectral line has its own values for /3 and 7, so we need 
to know the "effective" values of those line parameters in order to make a comparison between 
the analytic formula and the simulated results. For the purpose, we employ an autocorrelation 
method. By autocorrelating the synthetic spectrum for a GO-type star, we obtain an effective line 
profile of the spectrum (see AppendixOfor details). Since there is a strong degeneracy between the 
two stellar line parameters (3 and 7, we fix (3 based on a simple physical principle and estimate the 
Lorentzian dispersion 7 from the effective line profile. In principle, the intrinsic Gaussian dispersion 
/3o in each spectral line is determined by the effective temperature T e g of the star and the "the 
micro-turbulence" dispersion £ as 



po = f-^+e, (19) 

where k-g an d u are the Boltzmann constant and the mass of the atom (or molecule) in question, 



respectively (|Grayil2005l ). In addition, as we will show later, the simulated velocity anomalies also 
depend on the instrumental profile which we assume in generating mock spectra. Therefore, we 
here adopt an ad hoc assumption that the Gaussian width parameter (3 in Equation (|16p depends 
also on the width of the instrumental profile so that f3 is expressed as 



p = + a 2 p = \l 2J ^ +e+ /?i 2 p, (20) 

where /3rp is the dispersion of Gaussian which best-fits the representative instrumental profile 
adopted in observations. 



Sub stituting T pff = 600 K = 0.517 eV, p = 52 GeV/c 2 (the mass of an iron atom_j), f = 1-0 



km s 1 ( Coelho et al.ll2005l ). and /3rp = 3.6 km s 1 (the dispersion of the instrumental profile 



4 Iron lines are most numerous in the wavelength range used for the RV analysis with the iodine cell. 
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of Subaru/HDS we assumed in generating mock transit spectra, corresponding to the spectral 
resolution of R ~ 45000), we obtain /3q = 1.7 km s _1 and /3 = 4.0 km s _1 . Fitting the effective line 
profile with the Voigt function assuming /3q = 1.7 km s _1 , we obtain 7 = 0.9 km s _1 (Appendix 

ED- 

In each panel (different usini s ) of Figure [21 the analytic curves based on Equation (I16p are 
shown in the same colors as the simulated results indicated by symbols for each value of /. We 
adopt (/3,7)=(4.0 km s _1 , 0.9 km s _1 ). The other parameters (ui,U2,vsini s ,Q in Equation (|16p 
are fixed at the same values used to make mock transit spectra. 
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Fig. 2. — Simulated velocity anomalies due to the RM effect v.s. the analytic formula (Eq. [16]) as 
a function of the subplanet velocity v p in the case of a GO-type star: vsi~ni s = 2.5 km s _1 (top), 5.0 
km s _1 (middle), and 7.5 km s _1 (bottom). In each panel, symbols in black, red, blue, purple, and 
green indicate the simulated data points for / = 0.004, 0.008, 0.012, 0.016, and 0.020, respectively. 
Analytic curves based on Equation ()16|) are plotted with the solid lines in the same color as the 
symbols. 
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For all of the three different rotational velocities, the analytic curves are in very good agreement 
with the simulated results. As long as v sini s i s small, the velocity anomaly curves are nearly linear 
and can be approximated as Av ~ —fv p . As IWinn et al.l (|2005l ) and other authors have pointed 
out, however, they are significantly curved for larger values of vsini s . 



3.1.3. The Impact of Instrumental Profiles 

In the forward modeling method using the iodine cell, radial velocities v^y of a star are 
computed by modeling each observed spectrum iobs(A) with the following equation: 

J obs (A) = k [A(\)T (A(l - vrv/c))] * IP, (21) 

where ^4(A) and T(A) are th e transmission spectrum of the Iodin e cell and the template spectrum of 



the same star, respectively (jSato et al.ll2002l : iButler et al.lll996l ). Outside of a transit, the intrinsic 



stellar line profile of each spectrum is supposed to be the same as that of the template spectrum, 
and therefore urv is not affected by the instrumental profile which may be variable during an 
observation. During a transit, however, the intrinsic stellar line profile is distorted due to the 
partial occultation, and it is not clear if the IP affects the radial velocity urv including the RM 
velocity anomaly Av. Thus, in order to see if the instrumental profile used in the mock data analysis 
affects the RM velocity anomaly, we repeated the same mock data analysis described above, but 
with a different instrumental profile. In the simulation above, we have fixed the instrumental profile 
so that it corresponds to the spectral resolution of R = 45000 for Subaru/HDS. This time, we adopt 
the instrumental profile with the spectral resolution of R = 90000. In this case, the dispersion of 
Gaussian fits the instrumental profile is approximately /3rp = 1.9 km s . After generating many 
mock spectra for vsvai s = 2.5, 5.0, and 7.5 km s" 1 as in Section f3-l.lt we analyzed them with the 
usual RV routine to obtain the velocity anomaly Av. 

In order to quantify the differences between the two cases of the instrument profile, we compute 
the following statistics for each case of the instrumental profiles: 



J_ I Av sim (f, v p ) - At; ana (/, v p , (3, 7) 



\ 315 ^— ' I f ■ v sini s 

\ f,Vp,v sin i s 



(22) 



where Av s i m {f, v p ) is the simulated velocity anomaly for each point of the (/, v p ) grid, and Au ana (/, v p , f3, 7) 
is the value computed by Equation ([To]) as a function of /, v p , /3, and 7. When we compute 
Av a , n3i (f, v p , (3, 7), all the parameters other than f3 and 7 are fixed at exactly the same values as are 
adopted in the numerical simulation (such as £ and U\, 112)- The summation in the above expression 
is performed over the 315 (= 105 x 3) data points on the (/, v p ) grid with three different values of 
i>sinz s . The statistics, D({3,j), indicates the degree of agreement between the analytic formula and 
the simulated results in terms of the representative velocity anomaly fv sin i s . For each of the instru- 
mental profiles {R = 45000 and R = 90000), we compute D(/3, 7), for 2.0 km s" 1 < /3 < 5.0 km s" 1 
and 0.0 km s _1 < 7 < 3.0 km s~ . 
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Fig. 3. — The contour of -D(/3,7) for a GO- type star. The regions surrounded by the two blue 
and the two red curves are the best-fit regions of the analytic formula, where D(/3,7) < 0.005 for 
R = 45000 and R = 90000, respectively. The values of (/3, 7) estimated by the line analysis are 
shown by the blue and red crosses for R = 45000 and R = 90000, respectively. 

As a result of computing D((3,^/) in the grid, we find the lowest value of D(/3,7) is 

approximately 0.0045 for both of the two different instrumental profiles, representing very good 
agreement. For instance, the dispersion of the simulated velocity anomalies around the analytic 
formula is expected to be less than ~ 0.25 m s _1 in case of / ~ 0.01 and v sini s ~ 5.0 km s" 1 . 
This deviation is much smaller than the usual precision with which RVs can be measured. 

In Figure [3l we show contours of the goodness-of-fit statistic D(/3,7). The regions between 
the two blue curves and two red curves are where D(f3, 7) < 0.005 for R = 45000 and R = 90000, 
respectively. Since the Gaussian and Lorentzian dispersions f3 and 7 strongly degenerate, the 
values of f3 and 7 that fit the simulated results well are expected to be located in extended areas. 
Thus, we show the "best-fit regions" where the analytic formula agrees well with the simulated 
results. From Figure O it is obvious that the best-fit region is shifted toward smaller values of 
(/3,7) when we adopt the higher spectral resolution. This result implies that the instrumental 
profile in spectroscopic observations does affect the velocity anomaly due to the RM effect. We 
interpret the results as follows; in the forward-modeling fitting procedure, the instrumental profile 
(a finite spectral resolution) plays a similar role as the other physical broadening mechanisms for 
spectral lines. For reference, we show (/3,7)=(4.0 km s , 0.9 km s _1 ) by the blue cross, which 
are the intrinsic line parameters for R = 45000 estimated by the spectral line analysis in Appendix 
ICl Also, substituting /3rp = 1.9 km s" 1 into Equation (|20p for the case of R ~ 90000, we obtain 
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(3 = 2.5 km s~ x , which, along with 7 = 0.9 km s" 1 (the same value as used in the comparison for 
R = 45000), is shown in Figure [3] by the red cross. It should be emphasized that the values of (/3, 7) 
shown by the blue and red crosses are estimated in a way independent of the mock data simulation 
(Appendix [C]) but are consistent with the regions where the analytic formula best agrees with the 
simulated results, for the two different instrumental profiles. 

In summary, as a result of trying two different cases of the spectral resolution (R ~ 45000 and 
90000), we have shown that the RM velocity anomaly is actually affected by the specific choice 
of the instrumental profile. Since instrumental profiles are often approximated as Gaussian, we 
incorporate its impact on the line profile into the intrinsic Gaussian dispersion /3q by adding the 
instrumental broadening /3rp in quadrature, which justifies the treatment shown in Equation (|20p . 



3.1.4- Results for Other Spectral Types of Stars 

So far, we have considered a GO star. In order to make sure that our analytic formula is 
applicable for a variety of different stellar types, we consider stars with the effective temperatures 
of T e g = 6500 K (an F5 star) and 5000 K (a K0 star) using the theoretically synthesized spectra, 
just as we did for a GO star. 



F5 star In generating mock spectra for an F-ty pe star (T e g- = 6500 K), we broaden the synthetic 
spectrum of an F5 star (jCoelho et al,ll2005l ) assuming the rotational velocity of vsini s = 5.0, 
10, and 15 km s^ 1 . The ot her adopted parameters include u\ = 0.32 and = 0.36 for the 
limb darke ning param eters ( Claretll2004l )- and Q = 6.2 km s _1 for the macroturbulent velocity 
parameter (|Grayll2005h , Except for those parameters, we perform exactly the same simulation 
described in Section [3.1.11 The simulated results for v s'mi s = 10 km s" 1 are shown in the 
upper panel of Figure IU along with the analytic formula (Eq. [16]) assuming 7) = (4.0 
km s _1 , 0.9 km s" 1 ). These Gaussian and Lorentzian dispersions are also estimated by the 
combination of Equation ([20]) adopting T c g = 6500 K, and the spectral line analysis using 
the auto-correlation method. Again, the analytic curves well describe the behavior of the 
simulated results. They are also in good agreement with each other for vsmi s = 5.0 km s" 1 
and 15 km s _1 , although for brevity we do not show all of those results here. 
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Fig. 4. — The comparison between simulated results (symbols) and Equation (|16j) (solid curves) 
for an F5 star with v sini s = 10.0 km s _1 {upper) and for a K0 star with usini s = 3.0 km s" 1 
{lower), respectively. 



KO star We also make the mock transit spectra for a K0 star (T e g = 5000 K) and put them 
into the RV routine. Since most of the K-type dwarfs are relatively slow rotators, we adopt 
small rotational velocities: v sini s = 1.5, 3.0, and 4.5 km s _1 . We fix the limb-dar kening 
param eters and the macro turbulent velocity parameter to be u\ = 0.65, U2 = 0.14 (j Claret 
2004 ). and ( = 2.8 km s" 1 (jValenti et al.ll2005l ). We compare the simulated velocity anomalies 
with Equation (|16p . As an example we show the result for vsini s = 3.0 km s" 1 in the lower 
panel of Figured! for which we assume /3 = 3.9 km s" 1 and 7 = 1.1 km s , as estimated by 
analyzing the synthetic line profiles for GO and F5 stars. 
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All the numerical simulations show that our new analytic formula reproduces the simulated results 
within the current RV precisions. The agreement als o suggests that we c an vali d ate the previ- 
ously reported empirical relations for Subaru/HDS (i.e. jNarita et al, l l2009al lbl. [201o3 lbl:l Hirano et al 



20 111 ), which are based on the similar mock data simulations. 



3.2. Sensitivity of the Formula to Line Parameters 

In the previous subsection, we have shown that the analytic formula (Eq. [TB]) gives an accurate 
description of the simulated velocity anomalies during a transit as long as we adopt appropriate 
values of /3 and 7 for a given type of star. However, there are several practical issues that must be 
addressed before the analytic formula is applied to real data analysis. First, the macroturbulent 
velocity parameters that we assumed in both the application of the analytic formula and in the 
mock data simulations will not be known a priori for real stars. Second, we have used theoretically 
synthesized spectra to generate mock transit spectra, but actual intrinsic line profiles in observed 
spectra may differ from theoretical ones. Finally, although we have adopted representative instru- 
mental profiles in both of the analytic formula and simulation, the instrumental profile is generally 
dependent on many factors such as temperature variations, the position on CCD, etc, even if we 
adopt the same spectrograph setups for observations. Therefore, it is important to investigate the 
sensitivity of Equation (I16p to those parameters affecting the line profile (including the instrumen- 
tal profile). In this subsection, we perturb the values of those line parameters and examine the 
resulting changes to the outputs of the analytic formula. 



3.2.1. Dependence on Macroturbulence 
First we consider changes in the macroturbulence parameter. The macrot urbulence dispersio n 



£ for a GO-type star like HD 209458 is empirically estimated as C = 4.3 km s" 1 (IValenti et al.ll2005h . 
We here change it by ±30% with respect to that value (i.e. C = 3.0 and 5.6 km s _1 ) and plot the 
velocity anomaly curves for HD 209458 assuming the other line parameters as (/3,7) = (4.0 km 
s _1 ,0.9 km s" 1 ) and the stellar rotational velocity of v sini s = 4.5 km s . The top panel in Figure 
[5] shows the three cases of the macroturbulence dispersion: C = 4.3 km s" 1 (black), 3.0 km s _1 (red), 
and 5.6 km s" 1 (blue). 
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Fig. 5. — Variations of RM velocity anomaly curves for HD 209458 (top), XO-3 (middle), and TrES- 
1 (bottom) based on Equation (|16p . In each panel, we change the macroturbulence parameter £ 
by ±30% (red and blue curves) from the literature-based values (black curves). The residuals of 
the red and blue curves from the black curves are also shown at the bottom for each system. For 
the spin-orbit mis alignment angle A, we assumed A = 0° for HD 209458 and TrES-1, and A = 37.3° 
(Winn et a.l. 2009"). resnectivelv. 
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Inspection of Figure [5] shows that the RV difference between the curves is at most ~ 1 m s" 1 
(~ 1.5% of the representative velocity anomaly (R p /R s ) 2 v sini s for HD 209458), which is less than 
the RV precision from an observational point of view. Therefore, the choice of £ for this type of 
star does not significantly affect the results. 

We also show the comparisons for other two systems: XO-3 and TrES-1. For an F5 star XO-3, 
we adopt f3 = 4.0 km s _1 , 7 = 0.9 km s _1 , and vsini s = 18.5 km s _1 , and try three different cases 
of 4.3, 6.2, and 8.1 km s _1 (changed by ±30% around the central value). Likewise, we compute 
the RM velocity anomaly for TrES-1, assuming (3 = 3.9 km s _1 , 7 = 1.1 km s~ 1 and v smi s = 1.5 
km s _1 . We choose values for the macroturbulence dispersion of £ = 2.0, 2.8, and 3.6 km s _1 . The 
results are shown in the middle and bottom panels in Figure The RV differences of the colored 
curves from the black ones are at most «5m s _1 (~ 3.5% of (R p /R s ) 2 v smi s ) for XO-3 and ~0.1 
m s _1 (~ 0.4% of (Rp/Rg) 2 vsinis) for TrES-1, respectively. These values are less t han the RV 



preci sion that has b een achieved for ea ch system, which are approximately 8ms 1 ( Winn et al 



2009J) and 10 ms -1 (jNarita et al.l 120071 ). respectively. 



3.2.2. Dependence on (3 and 7 



Next, we focus on the Gaussian and Lorentzian components of line profiles. Fixing the values 
of macroturbulence parameter £, we change f3 and 7 in Equation (|16|) . and plot the RM curves 
for the three systems above. Figure [6] presents the three different cases of 7) for HD 209458. 
We change both of /3 and 7 by ±30% from the values estimated by the spectral line analysis 
(Appendix [C]) while we fix £ at 4.3 km s" 1 . According to Figured! the largest discrepancy of the 
two colored curves from the black one is < 2 m s" 1 (~ 3.1% of (R p /R s ) 2 v sinz s ). Although this is 
still comparable to or less than the RV precision for this type of stars, it is slightly larger than the 
difference in the curves for different £ (see the top panel of Figure [5]) . 
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Fig. 6. — The variation of RM velocity anomaly curve for HD 209458. We change the Gaussian 
and Lorentzian dispersion (/3, 7) in Equation ()16p by ±30% (red and blue curve) from the values 
based on the line analysis (black curve). The residuals of the red and blue curves from the black 
one are also shown at the bottom. 

We also compare the analytic RM curves for different values of (/3, 7) in case of XO-3 (F5-type 
star) and TrES-1 (KO-type star). As a result, the differences in the RM curves are at most < 3 m 
s^ 1 (~ 2.1% of (Rp/R s ) 2 v sin i s ) for XO-3 and < 0.5 m s" 1 for TrES-1 (~ 2.0% of (R p /R s ) 2 v smi s ). 
These deviations are less than the current RV precision for each type of star. For a rapidly rotating 
star such as XO-3, the line profile is mainly determined by the stellar rotation and macroturbulence, 
and thus the intrinsic line profiles (thermal and Lorentzian broadening) are likely to be less impor- 
tant. On the other hand, for slowly rotating late-type stars as TrES-1, the intrinsic line parameters 
(/3, 7) are important for describing the line shapes. However, since the RM velocity anomaly is well 
approximated as Av ~ — fv p for narrow line profiles (slow rotators) and the velocity amplitude 
due to the RM effect is comparably small, an inaccurate estimation for (/?, 7) or a variation of the 
instrumental profiles (less than ~ 30%) are not important for slowly rotating stars. 

In summary, the relative importance of the various line parameters depends on the stellar 
rotational velocity and the stellar type. For rapidly rotating stars (vsini s > 10 km s _1 ), the 
macroturbulence dispersion is more important than thermal and natural profiles, while the Gaussian 
and Lorentzian dispersions become more significant as for moderately rotating stars (3.0 km s _1 < 
vsinig < 10 km s _1 ). Neither of the effects is important for slowly rotating stars (vs'mi s < 3.0 km 
s-i). 
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3.3. Comparison with the Published "Calibrations" for Keck/HIRES 

In Section 13.11 we compared the analytic formula with simulated results based on the analysis 
routine for Subaru/HDS. We have confirmed that Equation (|16p well approximates the simulated 
results with deviations less than ~ 0.5% of the typical velocity anomaly scale (i.e. fvsini s ). It 
is very interesting to see if the new formula (Eq. |16j ) also agrees with the simulated results in 
previously published papers using some instrument other than Subaru/HDS. If they are shown to 
be consistent with each other, it suggests we no longer need to perform numerical case-by-case 
mock data simulations for RM analyses. If they do not agree, it may justify further investigation 
into which approach is more accurate. Here, we consider four systems for which the RM effect has 
been measured with Keck/HIRES: HAT-P-4, HAT-P-14, XO-3, and HAT-P-11. We compare Equa- 
tion (|16p with the published empirical relations based on mock data simulations for Keck/HIRES. 

As discussed in Section 13.1.31 the instrumental profile affects the velocity anomaly due to 
the RM effect. Therefore, in order to make a comparison between Equation (|16|) and the empirical 
formulae, we need to know the representative dispersion of the instrumental profile of Keck/HIRES. 
We here adopt /3ip = 2.8 km s _1 in terms of velocity (a representative dispersion for the spectral 
resolution of R = 65000). Since other parameters are intrinsic stellar properties, we adopt similar 
values of /3o, 7, and ( to the ones adopted for the comparison in Section [37X1 to draw the analytic 
RM curves. 

Table 2: Calibrations of the RM effect drawn from the literature, and choices of the stellar line 
parameters used in Figure [71 The stellar parameters 7, v sini<j) are expressed in km s . 



System Empirical Relation /3 7 £ v sin i s Reference 

HAT-P-4 Av = - /Wl.36 - 0.628 ( -^-Y] 3^2 (X9 4J 5.5 Winn et al. (2011 ) 

HAT-P-14 Av = -fv p [1.58 - 0.883 ] 3.3 0.9 6.5 8.4 Winn et al. (2011) 

XO-3 Av = -fv p [l.QU - 1.036 (^^j) 2 ] 3.3 0.9 6.2 18.5 Wnnetal (2009) 

HAT-P-11 At; = -fv p 3.2 1.1 2.5 1.5 Winn et al. (2010b) 



- 23 - 



in 

> 

< 



CO 

■o 

Ul 
CO 

DC 



> 

< 



1 ) 

eck Simulation 
Gaussian formula 
OTS formula 



Equation = 4 - 1 km s 



time [hour] 



50 



„ 



-50 - 



-100 



(c) Equation (1 6) (t = 6.2 km s 

Ke' 



eck Simulation 
Gaussian formula 
OTS formula 
HD3861 + rotation kernel 



XO-3 
v sin i s 
X = 37.3 deg 



18.5 km s 




40 




w 

> 

< 



w 
"to 

Z5 
"D 

(D 

cc 



E, 


30 r 


20 - 
10 - 


TO 





sidu 


-10 - 
-20 - 


CD 

DC 


-30 L 


-2 



time [hour] 



Equation ( 16 y? = 6 5 km s 



eck Simulation 
Gaussian formula 
OTS formula 




-0.5 
time [hour] 



CO 

E 



3 -4 



1 1 r 

(d) Equation (16 



i 1 r— 

2.5 km s ) 
eck Simulation 
Gaussian formula 




HAT-P-11 

v sin i s = 1.5 km s 

?i = 103 deg 




-1.5 -1 



-0.5 0.5 
time [hour] 



1 1.5 



Fig. 7. — The comparison between Equation (|16p (black) and the empirical relations in literatures 
(red) for (a) HAT-P-4, (b) HAT-P- 14, (c) XO-3, and (d ) HAT-P-11. For reference, analytic curves 
based on the Gaussian formula by iHirano et al.l (120101 ) and the OTS formula are plotted in blue 
and green, respectively. The residuals of the empirical curves and the previous analytic formulae 
from Equation (|16p are also shown at the bottom in each panel. 



Figure [7] shows the comparison between Equation (|16p and the empirical relations in literatures 
for the four systems. The empirical relations in the literatures are summarized in the second column 
of Table [2] and shown with red curves in Figure [71 Analytic curves (Eq. [16]) are shown in black. 
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The adopted stellar line parameters to draw the analytic curves are also summarized in Table [2j 
In Table [H /j, 7, C, and y sin i s are expressed in km s~ 1 . For reference, we plot the Gaussian 



formula by iHirano et al.l (|2010l ) (i.e. Eqn.|27|) in blue and the OTS formula in green, which is 
expressed as Av = —fv p /(l — /). In drawing the Gaussian formula, we employ the width j3 p of 
the intrinsic spectral line (in the absence of stellar rotation) approximated as a single Gaussian by 
Pp = \/ @ 2 + C 2 ] where /3 and £ are the values shown in Table EJ 

Regarding HAT-P-4, HAT-P-14, and HAT-P-11, for which the rotational velocities are rela- 
tively small, Equation (|16p and the numerical calibration formulas (red) agree within a few m s" 1 . 
For the HAT-P-11 system in particular, the two curves are almost indistinguishable. This is because 
Equation (|16p is well approximated by Av ~ — fv p for the case in which stellar lines are sufficiently 
narrow. On the other hand, the discrepancy between them for the case of XO-3 is more significant 
(~ 10 m s -1 ) in comparison with other systems although they are in much better agreement in 
comparison to the previous analytic curves based on the Gaussian and the OTS formulae. 

There are at least two possible explanations for the discrepancies between Equation f)16[> and 
the Keck simulated result. First, the line parameters (/3,7, C) adopted in Equation f)16|) to plot 
the analytic curves in Figure [7] may be significantly different from the true values for the system. 
Specifically, the macroturbulence dispersion £ is not well known for massive stars, and different 
values of £ lead to different results as we have shown in Section [3.2.11 In Section [3. 1[ on the other 
hand, we have compared Equation (|16p with simulated results adopting the same values of £ in 
Equation (I16p as the ones used to generate mock transit spectra. 

A second possible explanation for the discrepancy between the analytic formula and the Keck 
simulation is that the coupling between the stellar rotation and macroturbulence (which was ignored 
in the numerical calibrations) has a non-negligible impact on the results. In most cases, the mock 
data simulations for Keck/HIRES used real observed spectra to create mock spect ra during a 



transi t. For example, to create mock spectra during a transit of the XO-3 system, IWinn et al 



(|2009l ) began with the template spectrum of HD 3861, which is the same spectral type as XO-3 but 
has a smaller rotational velocity (F5V, v sini s = 2.8 km s _1 ). Then, they broadened that spectrum 
with a pure-rotational kernel to mimic the spectrum of XO-3. This is a good choice in the sense 
that the starting spectrum already involves the macroturbulence, which is not well known for those 
types of stars. However, as we have emphasized, the rotational and macroturbulence broadenings 
are coupled and, strictly speaking, cannot be applied sequentially. To be more precise, we should 
first deconvolve the spectrum with the macroturbulence before convolving with a kernel including 
the effects of both rotation and macroturbulence. 

In order to test the second hypothesis, we perform another mock data simulation for XO-3. 
We begin with the F5-type synthetic spectrum and broaden it with M(v) in which we assume 
«sini s = 2.8 km s" 1 and £ = 6.2 km s , so that the resulting spectrum mimics the template 
spectrum of HD 3861 (let us call it Ihd386i(^))- Then, we additionally broaden Thd386i(^) with 
the pure-rotational broadening kernel without macroturbulence assuming vsmi s = 18.5 km s _1 (we 
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cal l the resulting spec trum as 7xo3,Keck(^))- This process is expected to reproduce the steps taken 
by I Winn et al,l (|2009l ) in creating the mock spectrum of XO-3, although our starting spectrum is a 
theoretically synthesized one rather than that of HD 3861. By changing (f,v p ), we generate many 
mock transit spectra from 7xo3,Keck(A) and put the mock spectra into the RV analysis routine. 



As a result of fitting the output velocity anomaly Av as a function of / and 
following empirical relation: 



v p , we obtain the 



Av = -fv f 



1.60 - 0.900 



18.5 km s _1 



i crflt, crflt = 0.0073 km s" 



(23) 



where a&t is the dispersion of the residuals of the simulated Av from the best-fit curve. We plot 
this empirical relation in Figure [7] (c) in purple. Although we still see differences between the Keck 
simulation (red) and our new simulation (purple), they are reduced, and the new empirical relation 
follows more closely the variation seen in the published numerical calibration. This suggests that 
part of the discrepancy (and perhaps most of the discrepancy) between Equation (fTBT) and the Keck 
empirical formula is ascribed by the coupling between the stellar rotation and macroturbulence 
which was neglected in performing the numerical calibration. Note that the RV difference between 
the two simulations (red and purple) is comparable to the dispersion of our simulated results 
(fjflt = 7.3 m s _1 ). 



4. Impact of Differential Rotations of Stars 



One of the important applications of the new ana l ytic fo rmula is describing the stellar dif- 
ferential rotation in RM measurements. iGaudi <fe Wind (J2007J) previously pointed out the impact 
of differential rotations via RM measurements is negligible since they are beyond the precision of 
the RV measurements. However, their argument on the detectability is focused on relatively slow 
rotators. Since the RM effect is now measured for fairly rapid rotators (vsmi s > 10 km s _1 ), it is 
important to see if the impact of stellar differential rotations is still negligible based on the precise 
modeling of the RM effect described in this paper. In this section, we describe the impact of stellar 
differential rotations in the velocity anomaly curves. 



Following iReinersI (120031 ). we model the stellar angular velocity O, as a function of the latitude 
I on the stellar surface (I = at the stellar equator) as 



fi(i) = neq(l - a sin 2 /), 



(24) 



where f2 eq is the angular velocity at the equator. The coefficient a is the parameter describing 
the degree of differential rotation and approximately 0.2 for the case of the Sun. If the transiting 
planet is located at (x, y) on the stellar disk, its latitude I is estimated by 



■ , V ■ ■ 

sin I = — sin i s + 
it., 



x 2 + y 2 



R 2 



COS l s . 



(25) 
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When the host star is differentially rotating, we must replace the constant in Equations (jSJ) 
and (|B7p with Equation (|24p . In this case, the degeneracy between £1 and sini s can be solved; for 
rigidly rotating systems, the velocity anomaly Av depends solely on f2sini s , but with a differential 
rotation, the specific choice of the stellar inclination i s also affects the result. 
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Fig. 8. — The RM velocity anomaly curves for the XO-3 system with and without differential 
rotation. The blacks lines in each panel indicate the case when the host star is rigidly rotating 
while the other colors are for the cases that the star is differentially rotating. In the left panel, 
the stellar inclination is fixed at i s = 90° while a is changed from 0.0 up to 0.4. The RM curves 
for various stellar inclinations are shown in the right panel along with the case for rigid rotation. 
Different colors indicate different coefficients (a) and stellar inclinations (i s ) as shown in each panel. 
For each panel, RV residuals from the rigid rotation (the black line) are shown at the bottom. 



In order to evaluate the impact of differential rotations, we take XO-3 again as an example. 
Since the XO-3 system is reported to have a large spin-orbit misalignment and its host star has a 
large rotational velocity, we can expect a large impact of the differential rotation on RM velocity 
anomalies. Figure [8] indicates the RM anomaly curves for XO-3 with and without differential 
rotation. In the left panel, we change the differential rotation coefficient a, fixing the stellar 
inclination i s . We try the case of a = 0.0, 0.1, 0.2, 0.3, and 0.4. The right panel of Figure M 
indicates the variation of the RM velocity anomaly for the different cases of i s , in which we set 
i s = 90°, 60°, and 30° while a is fixed at 0.2. In this figure, the stellar spin velocity at the equator 
is fixed as v sini s = R s Q, cq smi s = 18.5 km s~ 1 . At the bottom of each panel, we show the residual 
of each curve from the black line (rigid rotation). 

It is notable that the residuals are comparable to or rather larger than the RV precision for 
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this system, which is reported to be ~ 8 m s _1 ( Winn et al.1 12009 ). The largest impact (~ 15 m 
s _1 ) on Av occurs for the case of large values of a and small values of i s , as is expected. Figure [8] 
also shows that the stellar inclination i s significantly affects the velocity anomaly. This result, in 
turn, suggests that we may be able to estimate the stellar inclination i s from RM analyses if we 
can model the differential rotation for each type of stars and fix the coefficient a. 



5. Discussion and Summary 



We have developed a new analytic formula (Eq. [TB]) to describe the velocity anomaly during 
a planetary transit. We here summarize the previous findings describing the RM velocity anomaly 
along with our new analytic formula: 

/ 



Av 



Av 



Av 



1 



(OTS) 

3/2 

fVp 



1 



+ 



(jfflrano et alJ 1201(1 ) 



/ 



Pl + Pj 2(/?2 + /3 p 2 ) 2 
exp(-2vr 2 /3 2 cr 2 - ^a)M(a)Q(a) sm(2Trav p )ada 



2tt 



{1 
M(a) - /0(er) cos{2irav p )j a 2 da 



(26) 
(27) 

(this work) 



where /3+ and /3 P in Equation (|27p indicate the best-fit Gaussian dispersions for the stellar template 
spectrum (including stell ar rotation), and for the intrinsic line broadening in the absence of stellar 
rotation, respectively (see lHirano et al.1 (|2010l ) for details). Major improvement in our new formula 
is that it precisely modeled various effects which can possibly affect stellar line profiles such as 
rotational broadening, macroturbulence, and the instrumental profile. 

Numerical simulations using the mock transit spectra and the RV analysis routine for Sub- 
aru/HDS have shown that our analytic formula agrees with the simulated results within ~ 0.5% of 
typical RM amplitudes (< 0.25 m s _1 in the case of a typical transit of a hot Jupiter). 

We have also compared Equation (|16p with the output of numerical calibrations of the RM 
effect that have been given in the literature, based on mock data simulations for Keck/HIRES. 
We have plotted the RM velocity anomaly curves for four existing transiting systems: HAT-P-4, 
HAT-P-14, XO-3, and HAT-P-11. Equation (|16p proved to be in good agreement with previously 
reported empirical relations, although some deviations were seen for the XO-3 system. We suggest 
that this disagreement comes from the coupling between rotational and macroturbulent broadening, 
and perhaps also to imperfectly known stellar line parameters. One possible means by which this 
issue could be clarified is to apply our new analytic formula to observed RVs and see if Equation 
(fTBT) gives a better fit to the data. We do not attempt to do so in the present paper, but we intend 
to reanalyze many of the existing data-sets in due course. 



As lHirano et al.1 (J201Cj) and others have pointed out, the relation between the position of the 
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planet and the observed RM velocity anomaly usually does not significantly change the estimation 
for the spin-orbit misalignment angle A since a spin-orbit misalignment is most likely to be shown by 
an asymmetry in the velocity anomaly curve during a transit. However, when the transit impact 
parameter b is small, the asymmetry in the velocity anomaly curve becomes less clear and the 
amplitude of a series of Av becomes more important. 

Analytic approaches to a phenomenon give us an insight into the behavior of the phenomenon, 
and sometimes reveals important aspects which numerical approaches are likely to overlook. We 
hope that our new formula will be useful in data analyses of the RM effect and make a contribution 
to the progress in this field. 
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In this appendix, we derive Equation (fT6j) based on Equations (pQ), (fTUl) . (fTTj) . and (fT2"j) . 

In order to calculate convolutions and cross-correlations, dealing with them in the Fourier 
domain, where they are expressed as products, can significantly facilitate the computation. First, 
we Fourier-transform Equations ([I]) and (llOj) as 



A. Derivation of the Analytic Formula fl!6|> 



-^"star {c ) 
J~ transit (c ) 



S(a)M(a), 

S(a)M(a) + fS(a)Q(a)e 



(Al) 
(A2) 



where the Fourier transform of an arbitrary function F(v ) is defined as 




dv 



(A3) 



so that M(a) and 0(a) are explicitly given by Equations (|17|) and f)18[) . respectively (see Appendix 
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IB"]) . Thus, the Fourier transform of the cross-correlation function C(x) becomes 

C(a) = -Fstar (cr) • transit (<7 ) 

= {§(a)Y M{a) \m(<j) - fQ(a)e- 2niav A . (A4) 
The Fourier transforms of the Gaussian and Lorentzian functions respectively are written as 

G(v) = -± r e-v 2 /? 2 =► G{a) = e~^\ (A5) 
P v 7r 

L W = --T^ =► L(a)=e- 2 ^H. (A6) 

When we assume the intrinsic line profile as S(v) = V(v;f3, 7) = G(u : /?) * L(t>;7) (the Voigt 
function), the Fourier transformation of S{v) becomes 

S(a) = e ^W-2*7M. ( A7 ) 
In this case, the inverse Fourier transformation for C(o) is expressed as 

/oo 
-00 

= / exp(-27r 2 /3 2 (T 2 - 4tt7H + 27r«7a?)M(<7) (m(ct) - /9((j)e~ 27rifT,, f ) da 
j —00 

= 2y exp(-27r 2 /3 2 cr 2 - 47T70-) |M(cj)| cos(2vrcra;)(i<T 

poo 

-2/ / exp(-27r 2 /3V - 4ir"/a)M (a)Q (a) cos {2ir a (v p - x)} da. (A8) 
■/ 



The derivative of Equation (|A8j) with respect to x is 
dC(x) 



/"°° r 1 2 

= -47ry exp(-27r 2 /3 2 cr 2 - 47T70-) |M(o-)| <t sin(27ro-x)(ia- 

poo 

-4vr/ / exp(-2vr 2 /3 2 cr 2 - 4vr7cr)M((j)e((T)<Tsin{27rcr(vp - x)}dcr. (A9) 
The maximum of the cross-correlation function is at x = Av, where 

exp(-2vr 2 /3 2 cr 2 - A^a) {M{<j)^ a $m(2it<7Av)d(j 

POO 

= -// exp(-27r 2 /3 2 a 2 - Airja) M (a)B (a)a sm{2Tra(v p - Av)}da. (A10) 
Jo 

Because of the exponential factors, the integrals in Equation ()A10p can be cut off at a ~ 1/2/3. 
Moreover, since the Gaussian width parameter /3 is of the order of several km s _1 and the velocity 
anomaly is typically less than one hundred m s" 1 , we obtain 

2iraAv < 0.1, (All) 
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and we can safely approximate as sin(27rc7Au) ~ 2naAv. Thus, using the identity for trigonometric 
functions: 

sm{2ircr(vp — Av )} = sm.(27rav p ) cos(2ivaAv) — cos(2Trav p ) sm(2waAv) 

Ri s'm(2irav p ) — 27tctAv cos(27T(TVp), (A12) 

we obtain following analytic formula for Av: 

/•OO 

/ exp(-27r 2 /3 2 (j 2 - 4n7a)M(a)@(a)sm(2irav p )adcr 

Av^-f Jo . (A13) 

J exp(-27r 2 /3V - A^a)M{a) |M((t) - f@{a) cos(2irav p )j a 2 da 



Note that contrary to previously reported empirical relations, Equation (|A13P implies that the 
velocity anomaly due to the RM effect is not proportional to the loss of flux /, but the second 
order of / slightly affects the results. 



B. The Fourier Transformation of the Rotational and Macroturbulence Kernel 

One of the great successes in dealing with the RM effect in the Fourier domain is that the 
rotational and macroturbulence broadening kernel M(v) becomes significantly simpler in the Fourier 
domain. Here, we Fourier-transform M(v) and derive Equation (I17p . 

Since the velocity component v appears only in the macroturbulence kernel 6(f), we can write 
the Fourier transformation of M(v) as 

/oo 
M(v)e- 2niav dv 
-oo 

/oo 
&(v - xQ sin i s )e- 2niav dv. (Bl) 
-oo 

The last expression above is expressed by a linear combination of the following form of an integral: 

00 1 f fv - xttsini s \ 2 ) , . . 

exp { - - 2mov \ dv, (B2) 



where k is a constant. The exponent in the above integrand reduces to 

v-xQ,smi s \ 2 . (v - xVL sin i s + mK 2 a) 2 2 . . . 

Ziriav = g \- irai(irK cri — 2xllsmi s ). (B3j 



K K 



Thus, using the Gaussian integral: 



-^e~ {v/K)2 dv = 1, (B4) 
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Equation (|B2p becomes 

exp < — ( ) — 2niav \ dv = exp{iria(iriK 2 a — 2x0, sin i s )}- (B5) 



00 1 / v — xQsmie^ 



Replacing k with £cos# or £sin#, the second line in Equation (IBip is expressed as 

/oo — 2nicrx£l sin i s 
6(u - xfisini^e- 2 ™ 1 '^ = | e -Kcose)V 2 + e -«sin0)V|_ (Bg) 
-00 ^ 

Therefore, we obtain the following expression for M(cr): 

^ _ rr 1 - ui(l - cosfl) - it 2 (1 - cosfl) 2 

J Jen 



' entire disk 



7r(l - 14i/3 - U2/6) 



„— 2niaxQ sini s „ „ „ „ 

x £ l| e -«eos6?)V 2 + e -«sine)V 2 |^y _ (B?) 

Since the imaginary part of the integrand in Equation ()B7p is an odd function with regard to x, 
it vanishes after integrated over x. We then change variables as x = rcoscj) and y = rsin^, which 
leads to 

rR 3 r 2n 1 _ Ul Q - COS 0) - u 2 {l ~ COS 6) 2 



M(a) 



Jo 



7r(l - 14i/3 - IX2/6) 



-« cos 9) 2 cr 2 , -(7rCsin6l) 2 o- 2 



x cos(27r<7rf2sini s cos 0)- 5 — , (B8) 



where 



cos6 = Jl-j^, s[n6 = Y- (B9) 



Using the following formula 

»2tt 



/"Z7T 

/ cos(27ro"rO sini s cos 4>)d(j) = 2TrJo(2TiarQ sini s ), (BIO) 
Jo 



we obtain 

Af(o-) 
or equivalently, 

M(cr) 



Rs l-i4i(l-cos6')-i42(l-cos6') 2 



1 — 14i/3 — 1*2/6 

x {e -«eose) 2 CT 2 + e-W^T^jj^arnsinis) 7 ^-, (Bll) 



1 1 - 1x1(1 - ^/T^) - u 2 (i - VT^tf) 2 
l — 1x1/3 — 142/6 

x (e-^C'^C 1 -* 2 ) + e — 2 C 2 - 2 * 2 } J (2irav smi s t)tdt. (B12) 
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C. Estimation of the Effective Line Profile 

Here, we describe how to estimate the "effective" line profile from a spectrum. First, we 
model the auto-correlation function of a spectral assuming a simple analytic form of it, and then 
compute the actual auto-correlation function for the synthetic spectra that we use in the mock data 
simulation. 



We assume that the synthetic spectrum (jCoelho et al.ll2005l ) is expressed by a summation of 
Voigt functions (defined by Eq. [13]) as 

# of lines 

T(A)= Yl kV(\-di;f3 ini ), (CI) 

i 

where hi, di are the equivalent width and the central wavelength of the spectral line i. The Gaussian 
and Lorentzian dispersions of the line i are also noted by and 7 j, respectively. For simplicity, we 
subtract the continuum (normalized to unity) of the spectrum and flip the absorption lines so that 
the flux takes positive values. The auto-correlation function A(x) of the spectrum is defined as 

"OO 



/oo 
T(X)T(X-x)d\ 
-00 

/oq # of lines # of lines 
]T ]T kikjViX-diifcriViX-dj-xiPj^dX. (C2) 
i j 

Since the convolution between two Voigt functions yields another Voigt function with a larger 
dispersion, the integral above reduces to 

# of lines # of lines 



A ( x ) E hkjVix-di + djiJpf + pi^i + jj) 



* 3 

# of lines # of lines 



E kfV(x;V2f3 i ,2 li )+ £ £ kikjV(x - \d { - dA; Jff + /3|, 7 , + 7j ) 

# of lines 



+ E Y^kikjVix+ldi-dA^ff + pl^ + jj). (C3) 

i i>i 

The first term in the last expression indicates the sum of all the absorption lines broadened by 
factors of y/2 for Gaussian and 2 for Lorentzian profiles. Here, we use this term to calculate the 
"effective" line profile of all the spectral lines; we define the effective Gaussian and Lorentzian 
dispersions j3 and 7 as 

# of lines # of lines 

£ ^y(x;v / 2^,27)= Yl k?V(x;V2Pi,2 7i ). (C4) 
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Although a summation of Voigt functions is not necessarily expressed by a Voigt function, this 
treatment above is justified as long as all the values of /?, and ji do not differ significantly from line 
to line. 

The second and third terms in Equation (IC3j) mean the correlations between the lines separated 
by \di — dj\. It is plausible to assume that this separation \di — dj\ is distributed randomly in a 
spectrum. In this case, the second and third terms become the summations of a bunch of Voigt 
functions centered at randomly placed x. Assuming the number of spectral lines is large enough 
(i > 1000), the second and third terms altogether becomes almost constant as a function of x Q 
Therefore, we can approximate Equation (|C3P around x = as 

# of lines 

A(x) » W{x; V2fi, 2 7 ) + C, k = ^ kf, (C5) 

i 

where C is the constant originated from the second and third terms in Equation (]C3p . 

Next, we try to compute the auto-correlation function A(x) for the actual synthetic spectra 
we use to generate the mock transit spectra. We take the GO-type spectrum described in Section 
13.1-H subtract the continuum, and flip the lines. Then, we compute the auto-correlation function 
A(x) defined as Equation (|C2|) . Note that we cut the spectrum so that it ranges from ~ 5000 A to 
~ 6000 A, where precise radial velocity analyses are usually performed. 
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Fig. 9. — The auto-correlation function A(x) of the synthetic spectrum for a GO-type star, defined 
by Equation (|C2|) . The wavelength range used in the calculation is 5000 A < A < 6000 A. 

Figure [9] shows the resultant auto-correlation function A{x) of the synthetic spectrum for a G0- 




In reality, since we deal with a limited wavelength range of the spectrum, the auto-correlation function A(x) 
slowly decreases as a function of 
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type star used in the mock data simulation around x = 0. As is expected from the above theoretical 
prediction, it has a strong Voigt-shaped peak at x = with an offset due to the correlations with 
neighboring lines. We fit this auto-correlation function assuming a form of Equation (|C5p , Since 
the Gaussian and Lorentzian dispersion j3 and 7 are strongly correlated, we fix /3 as j3 = /3q = 1.7 
km s , which is calculated by Equation (|19p . and let only k, 7, and C be free. As a result, we 
find the best-fit value of 7 as 7 ~ 0.94 km s _1 . In this paper, we use these values in comparing the 
analytic formula with the simulated results for a GO-type star. 
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